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ABSTRACT 

This is the first of a series of papers investigating the oscillation properties of relativistic, 
non-selfgravitating tori orbiting around a black hole. In this initial paper we consider the ax- 
isymmetric oscillation modes of a torus constructed in a Schwarzschild spacetime. To simplify 
the treatment and make it as analytical as possible, we build our tori with vertically integrated 
and vertically averaged quantities, thus transforming the eigenvalue problem into a set of cou- 
pled ordinary differential equations. The tori are also modeled with a number of different 
non-Keplerian distributions of specific angular momentum and we discuss how the oscillation 
properties change when different distributions of angular momentum are considered. Our in- 
vestigation progresses by steps. We first consider a local analysis in Newtonian gravity and 
determine the properties of acoustic wave propagation within these objects, as well as the 
relations between acoustic and epicyclic oscillations. Next, we extend the local analysis to 
a general relativistic framework. Finally, we perform a global analysis and determine both 
the eigenfunctions and the eigenfrequencies of the axisymmetric oscillations corresponding 
to the p modes of relativistic tori. These behave as sound waves globally trapped in the torus 
and possess eigenfrequencies appearing in the simple sequence 2:3:4:. . ., independently of 
the distribution of angular momentum considered. The properties of the modes investigated 
here are in good agreement with those observed in recent numerical simulations and could 
have a number of different applications. In X-ray binary systems containing a black hole can- 
didate, for instance, p-mode oscillations could be used to explain the harmonic relations in the 
high frequency quasi-periodic oscillations observed. In systems comprising a massive torus 
orbiting a black hole, on the other hand, p-mode oscillations could be used to explain the 
development or the suppression of the runaway instability. 
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1 INTRODUCTION 

Waves and normal-mode oscillations in stars have long since at- 
tracted attention and a vast literature, investigating them both in 
Newtonian (see, for instance. Cox 1980, or Unno et al., 1989) and 
in general relativistic regimes (see Stergioulas 1998, for a review), 
is now available. On the other hand, waves and normal-mode os- 
cillations in geometrically thin discs around compact objects have 
been studied much less, both within Newtonian gravity (see Kato 
2001 for a review) and within a relativistic framework (Okazaki 
et al., 1987; Perez et al., 1997; Silbergleit et al., 2001; Kato 2001, 
Rodriguez et al. 2002). The literature is even more scarce when 
one considers geometrically thick discs, which so far have have 
been investigated mostly in connection with their stability proper- 
ties both in Newtonian gravity (Papaloizou and Pringle 1984, 1985; 
Blaes 1985) and in general relativity (Kojima, 1986). An explana- 
tion for why discoseismology has not yet reached the development 



and the level of sophistication which is now possible in asteroseis- 
mology is due, at least in part, to the fact that only recently ac- 
cretion discs have been recognized as fundamental astrophysical 
objects, present in various forms at all scales. Nowdays, however, 
periodic and quasi-periodic variations are currently observed in dif- 
ferent classes of astrophysical objects containing accretion discs. 
Although many different models have been proposed for the inter- 
pretation of the rich phenomenology associated with these quasi- 
periodic oscillations (QPOs), there seems not to be yet a widely 
accepted mechanism for most of the observed sources (see van der 
Klis, 2000 for a review). Clearly, a systematic investigation of the 
oscillation properties of discs could shed some light on this. 

As in stars, oscillation modes in discs are, in general, the 
consequence of restoring forces responding to perturbations and 
these offer a way for classifying oscillations. In accretion discs, in 
particular, a first restoring force is given by the centrifugal force, 
which is responsible for the appearance of the so called inertial 
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waves, tightly related to the orbital motion of the disc and hence to 
epicyclic oscillations. A second restoring force is offered by pres- 
sure gradients and the oscillations that are produced in this way are 
related to p modes and have close connections with the propaga- 
tion of sound waves in the perturbed fluid. A third restoring force 
is the gravitational field in the direction orthogonal to the orbital 
plane. If a portion of the disc is perturbed in the vertical direction, 
in fact, the vertical component of the gravitational field will pro- 
duce a harmonic oscillation across the equatorial plane with oscil- 
lation frequency equal to the orbital frequency. These oscillations 
are related to corrugation waves (see Kato, Fukue and Mineshige, 
1998, for an overview on disc oscillations). 

Two complementary approaches have been followed tradition- 
ally for studying the perturbations of equilibrium models of geo- 
metrically thin accretion discs. The first one is a local approach, 
and it has been used extensively to investigate the propagation of 
waves with the inclusion of many contributing physical effects, 
such as buoyancy, stratified atmospheres and magnetic fields (see, 
among the others, Lubow and Ogilvie, 1998). Being local, these ap- 
proaches derive dispersion relations in which the frequency of the 
perturbation is expressed as a function of the spatial position within 
the object and as the linear combination of different contributions, 
each related to the different physical effect (Kato, 2001). The sec- 
ond approach is a global one and it is based on the fact that, under 
suitable conditions, Eulerian perturbations to all physical quantities 
can be expressed as an eigenvalue problem, the eigenvalue being 
the frequency of oscillation. In practice, and for the simplest case, 
this amounts in solving a second order partial differential equa- 
tion once appropriate boundary conditions are provided (see Ipser 
and Lindblom, 1992; Silbergleit et al., 2001; Nowak and Wagoner, 
1991). 

This paper is devoted to both a local and a global perturbative 
analysis of axysimmetric modes of oscillation of relativistic tori in 
the Cowling approximation (i.e. in an approximation in which the 
perturbations of the spacetime are neglected; Cowling 1941). Fluid 
tori differ from geometrically thin discs mostly in having equilib- 
rium configurations in which the orbital motion is intrinsically non- 
Keplerian and in which pressure gradients play an important role, 
giving rise to an extended vertical structure. As a result, a con- 
sistent investigation of these configurations which would account 
for the coupling of the oscillations in the radial and vertical direc- 
tions requires necessarily a two-dimensional treatment involving a 
set of partial differential equations. While this problem is solvable 
with presently available techniques (this has indeed already been 
done for relativistic rotating stars; Yoshida and Eriguchi, 1997), 
it would require a massive use of expensive numerical calculations 
leaving little room to a physical interpretation. Furthermore, in con- 
trast with the case of relativistic stars, the eigenvalue problem for 
relativistic tori is in great part unsolved and resorting to a fully nu- 
merical solution at this stage would prevent one from appreciating 
most of the basic physics behind the oscillation modes of these ob- 
jects. For this reason, we will here concentrate on a simpler model 
for the torus which can be handled in great part analytically, and 
postpone the use of the two-dimensional fully numerical analysis 
to a subsequent work. 

The main simplification in the models discussed here is that 
the vertical structure of the tori is accounted for by integrating the 
relevant quantities along the direction perpendicular to the equato- 
rial plane. Doing so removes one spatial dimension from the prob- 
lem, which can then be solved integrating simple ordinary differ- 
ential equations. The local approach, in particular, will allow us 
to derive the local dispersion relation obeyed by the oscillations in 



relativistic non-Keplerian discs, while the global approach will pro- 
vide us with the eigenfunctions and eigenfrequencies of the system. 

There are several motivations behind this study. Firstly, we 
want to extend the relativistic discoseismology analysis carried so 
far for thin discs to systems having a non negligible contribution 
coming from pressure gradients. Secondly, we intend to interpret 
and clarify some of the numerical results found in the time evolu- 
tion of "toroidal neutron stars", i.e. compact and massive tori orbit- 
ing a Schwarzschild black hole (Zanotti, Rezzolla & Font, 2003). 
Thirdly, we want to assess the possible connections between the os- 
cillation modes of relativistic tori and the rich X-ray phenomenol- 
ogy observed in QPOs. Finally, we want to investigate the possibil- 
ity that the axisymmetric oscillations of thick discs could provide 
a criterion for the development or the suppression of the runaway 
instability (Abramowicz et al., 1983). 

The plan of the paper is as follows: in Section |2| we briefly 
review the local analysis of oscillation modes in vertically inte- 
grated Newtonian tori. Section |3| on the other hand, introduces 
the basic assumptions and equations employed in the definition of 
our general relativistic, vertically integrated torus. These equations 
will then be used to study axisymmetric oscillations both locally, 
in Section|4| and globally, in Section|5| We will first consider con- 
figurations with constant distributions of specific angular momen- 
tum, and subsequently distributions of specific angular momentum 
that are either linear or power-laws in the cylindrical radial coordi- 
nate. Finally, Section[8|contains our conclusions and the prospects 
for further investigations. Hereafter Greek indices are taken to run 
from to 3 and Latin indices from 1 to 3; unless stated differently, 
we will use units in which G — c — Mq = 1. 



2 NEWTONIAN TORI: A LOCAL ANALYSIS 

We here briefly present a local analysis of the oscillation modes of 
vertically integrated tori within a Newtonian framework. Some of 
the results presented in this Section have already been discussed in 
the literature, but will serve as a useful reference for the general 
relativistic treatment presented in Section |4| and will help in the 
physical interpretation of the relativistic results. 



2.1 Assumptions and Equations 

Consider an extended perfect fluid configuration orbiting a central 
object which is the only source of the gravitational potential (i.e. the 
orbiting fluid is non-self gravitating). Introduce now a cylindrical 
coordinate system, (uj, <f), z) whose origin is at the center of the 
central object and whose 2;-axis is oriented along the direction of 
the orbital angular momentum vector. 

A simplified description of this system can be obtained by re- 
moving the dependence of the various physical quantities on the 
vertical coordinate z. Mathematically, this is done by integrating 
the relevant physical quantities along the vertical direction. Phys- 
ically, this corresponds to collapsing the vertical structure of the 
torus onto the equatorial plane, but is quite different from just con- 
sidering an equatorial slice of the vertically extended torus. 

Using this approach, it is possible to define the vertically inte- 
grated pressure P 

P{zo) = f pdz, (1) 
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and the vertically integrated rest-mass density E 

T,{zu) = / pdz , 



(2) 



where H = H{vj) is the local "thickness" of the torus. An equa- 
tion of state (EOS) for the fluid in the torus needs to be specified 
and we find it convenient to use here a simple barotropic EOS of 
polytropic type, i.e. 



p — kp'' 



(3) 



where k and 7 = dlnp/dlnp are the polytropic constant and the 
adiabatic index, respectively. The vertically integrated quantities in 
Q and J2j need also to be related through an "effective" equation 
of state and this can be done if we define an "effective" adiabatic 
index T = d In P/d In E, so that 



P = ICE' 



(4) 



with the constants IC and F playing the role of the polytropic con- 
stant and of the adiabatic index, respectively. Note that while equa- 
tion J4} mimics a polytropic equation of state, it does not represent a 
vertically integrated polytropic equation of state (unless, of course, 
the adiabatic index is equal to one). More importantly, the adia- 
batic index F is not constant but depends both on uj and z. This 
complication, however, can be removed if one assumes that both 
the pressure and the rest-mass density have a weak dependence on 
height, so that they can be accurately expressed in terms of their 
values at the equator, i.e. 



: PO = p(ro, Z = 0) , 



p = p(rJ7, z) ^ pq = p{-aj, 2: = 0) 



(5) 



, - . (6) 

Using this assumption, P ~ 2Hpo and E ~ 2Hpo, so that 
F « dlnpo/dlnpo = 7 and equation J4j can effectively be writ 
ten as 



P = fcE^ 



(7) 



In all of the calculations reported here we have used F — 4/3, but 
the results do not change qualitatively when different polytropic 
indices are used. 

Being in circular non-Keplerian motion with angular velocity 
Q, the torus will have velocity components only in the vj and in the 
^-directions, with vertical averages given by 

f-H 



1 

2H 



v^dz 



and 



1 /-^ 



V dz . 



(8) 



(9) 



With these assumptions, the dynamics of the torus is fully de- 
scribed by a set of equations comprising the Euler equations in the 
TZ! and (/^-directions (Shu, 1992) 

W 1 

dtU + Ud^U + —d^U - = -_a„p_9^* ,(10) 



w . 



zo E 

uw 



dtW + Ud^W + —d4,W + 

-i-c»^p--a^*, (11) 

the equation for the conservation of mass 

atE + a„(E(7) + -Et/ + -a^(EW) =0, (12) 



and the Poisson equation for a vertically integrated gravitational 
potential ^ 



V * = 47r5 



(13) 



where S is the vertically integrated rest-mass density of the central 
object, only source of the gravitational potential. 

Stationary and axisymmetric models for the tori can be built 
after setting to zero the terms in jlO> - J12t involving derivatives 
with respect to the t and (^i-coordinates. The solutions obtained 
in this way represent the background equilibrium solutions over 
which harmonic Eulerian perturbations of the type 



I 5U \ 



(14) 



can be introduced. Note that Eulerian perturbations are indicated 
with a "5" to distinguish them from the corresponding Lagrangian 
perturbations, that we will instead indicate with a "A". In expres- 
sion J14t . a is the frequency of the perturbation (in general a com- 
plex number), fctc is the wavenumber in the ro-direction (hereafter 
simply k), and 5q = SP/E has been introduced to describe the 
perturbation in the pressure. 

We have here neglected the perturbations in the gravitational 
potential and therefore set S'if = 0. This is referred to as the Cowl- 
ing approximation (Cowling 1941) and for a fluid configuration 
which is non-selfgravitating (i.e. one whose gravitational potential 
is such that ^fluid + 5'1'fluid ~ 0), the Cowling approximation is 
actually an exact description of the pulsations (Ipser and Lindblom, 
1992). 

Note that the harmonic spatial dependence expressed in <14> is 
the signature of the local approach and is valid as long as the wave- 
length of the perturbations considered is smaller than the length- 
scale of the radial variations in the equilibrium configuration (this 
is basically the condition for the WKB approximation) or, stated 
differently, that 



k > 



(15) 



1 dP 1 dP 1 
P dzu P dzu L 
Introducing the perturbations <14t in the equilibrium model 
given by equations <10t - <13> and retaining only the first-order 
terms, we derive the following perturbation equations (Shu, 1992) 



iaSU + 2mW = ikSq , 



{2Q + wd.^Q)iSU + adW = , 



iaSq — ikc^SU = , 



(16) 



(17) 



(18) 



where = dP/dT, — FP/E is the local sound speed. We can 
now cast equations <16> - <18> into a simple matrix form as 



/a 20, -k \ 

2 



2n 







\ kc 







/ iSU \ 



5W 



(19) 



where is the epicyclic frequency in the radial direction^ and is 
^ In principle, an epicyclic oscillation takes places also for motions away 
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defined as 



azu 



(20) 



Note that the radial epicychc frequency is equal to the orbital fre- 
quency for Keplerian orbital motion, i.e. — for Q, ~ zu~^^^, 
and is zero for motions with constant specific angular momentum, 



i.e. for I 



const. 



Being a homogeneous linear system, a non-trivial solution to 
the set of equations <19t exists if the determinant of the coefficients 
matrix is equal to zero, which then provides the dispersion relation 



2 _ ^2 . , 2 2 



(21) 



This approximate form of the dispersion relation was first ap- 
plied to waves in accretion discs by Okazaki et al. (1987) and then 
reconsidered by several authors in more general situations (see 
Nowak and Wagoner, 1992; Ipser, 1994; Silbergleit et al., 2001). 
The two terms in the dispersion relation <21> are most easily inter- 
preted when considered separately. To this scope, consider a torus 
composed of collisionless particles and with specific angular mo- 
mentum increasing outwards. A fluid element which is infinitesi- 
mally displaced from its equilibrium orbit but conserves its angular 
momentum unchanged, will start oscillating in the radial direction 
due to a restoring centrifugal force. These oscillations are called 
inertial oscillations and their frequency is the radial epicyclic fre- 
quency Hir{'uj). In compressible fluids, on the other hand, a restor- 
ing force due to pressure gradients is also present and is responsible 
for acoustic oscillations with frequency kcs- Both of these terms 
contribute to the right-hand-side of the dispersion relation J2H and 
are collectively referred to as " inertial-acoustic waves" . Following 
a standard convention (Kato, 2001), we will here identify the high 
frequency modes (cr^ > ) with the p modes, or inertial-acoustic 
modes (Kato and Fukue, 1980) of the vertically integrated torus. 



3 RELATIVISTIC TORI: ASSUMPTIONS AND 
EQUATIONS 

As for the Newtonian case, we will here assume that the torus does 
not contribute to the spacetime metric, which we will take to be 
that external to a Schwarzschild black hole. This is the simplest, 
yet non-trivial metric to consider and, more importantly, it allows 
for a direct comparison with the numerical calculations of Zanotti 
et al. (2003) (The extension of this work to the Kerr metric will be 
presented in a future work; Rezzolla and Yoshida, 2003). Further- 
more, since we are really interested in the portion of the spacetime 
in the vicinity of the equatorial plane (i.e. for values of the spherical 
angular coordinate |S—7r/2| <C 1), we will write the Schwarzschild 
metric in cylindrical coordinates and retain the zeroth-order terms 
in the ratio {z/r) (Novikov and Thorne, 1973). In this case, the line 
element takes the form 



ds 



where 



„2i.(^) ,,2 



dt + e 



2A(tn)j 2 . , 2 I 2 



2^(i^) 



= 1 



2M 



-2A(tn) 



(22) 



(23) 



from the background orbital plane. The corresponding frequencies are the 
vertical epicyclic frequencies but these will not be considered here. 



The basic equations to be solved in order to construct equilib- 
rium models are the continuity equation for the rest-mass density p 



V<,(/9^t") = 0, 

and the conservation of energy-momentum 

V^T"'' = , 



(24) 



(25) 



where the symbol V refers to a covariant derivative with respect to 
the metric In equation |25j, T"'^ = {e + p)u°'u'^ + pg"*^ are 
the components of the stress-energy tensor of a fluid with isotropic 
pressure p and total energy density e. Since we are dealing with 
a curved background spacetime, but we want a close comparison 
with Newtonian expressions, it is useful to introduce an orthonor- 
mal tetrad carried by the local static observers and defined by the 
one-form basis 

a;* — dt , u>^ = e^dvj , ui^ = dz , uj"^ = Tudtf) . (26) 

In this frame, the components of the fluid four- velocity are denoted 
by u'^ (with fi = t,zu, z, 0), and the 3-velocity components are 
then given by 



V =^ = — , t = zu,z,(t). (27) 

In analogy with the Newtonian case presented in Section |2| 
we use the vertically integrated quantities defined in ™d 
assume that they obey an effective polytropic equation of state of 
the type J4}- Enforcing the conditions of hydrostatic equilibrium 
and of axisymmetry simplifies the hydrodynamical equations con- 
siderably, reducing them to Bernoulli-type equations (Kozlowski et 
al., 1978) 



e+p 



-Viln(Mt) + 



(28) 



where £ = —u^/ut is the specific angular momentum (i.e. the an- 
gular momentum per unit energy). Simple algebraic manipulations 
show that the only relevant component of equations <28t . i.e. the 
radial one, can be written explicitly as 

d^p _ _ e^^d^v - n^Tu ^29-) 
e + p e^"^ — Q?-co'^ 

which represents a generic condition for the existence of hy- 
drostatic equilibrium of a fluid configuration orbiting around a 
Schwarzschild black hole. Note also that a number of cancellations 
have removed any dependence on the spatial derivative of the an- 
gular velocity from the right-hand-side of equation <29> . 

Being interested in vertically integrated expressions, we need 
to integrate both sides of equation <29> . Interestingly, in the space- 
time <22> the right-hand-side of <29> is a function of the cylindri- 
cal radial coordinate only, and hence its vertical integration simply 
yields a new multiplicative factor 2H. The left-hand- side, on the 
other hand, cannot be managed unless a simplifying assumption is 
made and this amounts to consider that the total enthalpy (e + p) 
is only a weak function of the z-coordinate so that, effectively, it is 
possible to substitute it with its vertically averaged value (e + p), 
where 

1 

e{uj,z)+p{Tii,z) ^ e{uj)+p{uj) = — J {e+p)dz. (30) 

With this assumption, the vertically integrated equation <29> can be 
written as 

1 dP vj 



E + P dvj 



M 2 ^ 

— il ZU 



(31) 
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where E is the vertically integrated energy density 



E = 



edz = 



7 - 



(32) 



We next introduce perturbations in the velocity and pressure with a 
harmonic time dependence of the type 



\ SQ J 



(33) 



where SQ = 5P/ (E+P), and where we have defined the averaged 
velocity perturbations as 

SV"^ = -- / Sv'^dz , SV''' = -- Sv'f'dz . (34) 

After perturbing the hydrodynamical equations <24t - <25> and 
taking the vertically integration of the resulting equations, we de- 
rive the following set of equations 



E + P dvj 



(35) 



dvj zxj dvj J 

iaAne'^^vjSQ = , (36) 



E + P 
TP 



2^+ ' 



dzu 



A 



1 dA 1 dP 

+ 



dzo zu 2A dzu FP dzu 



e oV 



. 



(37) 

where A = l/(ut)^. Equations <35> - <37t represent the zu and ((>- 
components of the perturbed relativistic Euler equations as well as 
the perturbed continuity equation. 



4 PERTURBATIONS OF RELATIVISTIC TORI: A 
LOCAL ANALYSIS 

A local analysis of the perturbed hydrodynamical equations <35> - 
<37t can be performed after introducing a harmonic radial depen- 
dence of the type 5Q, SV^ , SV^ ~ exp (ikzu), where k is the ra- 
dial wavenumber and, again, we will assume that X = 2n /k <ti L. 
Simple considerations allow now to remove some of the terms in 
the continuity equation <37> . thus simplifying it further. In particu- 
lar, it is easy to realize that the third and fourth terms in equation 
<37t are very small as compared to the first and second ones, to 
which they are similar in nature. More specifically, both the sec- 
ond and the fourth terms involve the radial velocity perturbation 
but with coefficients proportional to k and 1/L, respectively. As- 
suming that the condition <14> is satisfied, it is then reasonable to 
drop the fourth term of equation <37> . Similarly, the ratio between 
the third and the first term can be approximated as 



(5V^ rPzune-" 
~5Q {E + P)A 



■ O 



5V^ 
~5Q 



X £>(! typical velocity n , (38) 



where we have approximated the first term as ~ 5P/{VP) ~ 
5Q/cl and the third one as ~ {vjQ,)6V^ . While {5V^ /5Q) is 
the ratio of two perturbed quantities and of order unity, the typical 
velocity of is of the order of the orbital velocity at the marginally 
stable circular orbit (i.e. ~ l/\/6) so that, overall, the ratio in <38> 
is much smaller than unity and the third term in equation <37> can 
therefore also be dropped. Finally, it should be noted that the third 
term on the left hand side of equation <36t can also be discarded, 
as it would introduce a purely imaginary part in the dispersion rela- 
tion <42t . This term, which is due a time derivative of the pressure 
in the Euler equations, has a purely relativistic origin and has been 
found to provide negligible contributions in the numerical solution 
of the system <35l-<37l. 

We can now introduce quantities that are more directly related 
to the ones used in the Newtonian example of Section|2| namely 



5W = 5V^ , 



(39) 



so that the linearized perturbation equations in the unknowns 5U, 
SW and SQ can be written as a homogeneous linear system 



/ SU \ 



M{k,a) 



SW 



\ 



(40) 



where the matrix of coefficients Ai depends both on k and a. Im- 
posing the determinant of A4 to be zero, we obtain the dispersion 
relation 



2n ( 2n + zu^ 

dzu 



dP 



E + P dzu 



+ e^ 



„ dv 

2zun— 

dzu 

-^"zu'^n^) e 



1 + 



E + P 



A non-trivial solution of <4U is given by 

^2 ^^2 + fe2(-^)(l_e-2- 



= . 
(41) 

(42) 



where now = dP/dE is the square of the relativistic sound ve- 
locity in the vertically integrated disc. The dispersion relation j42t 
represents the relativistic generalization of <21> in which we have 
defined the relativistic radial epicyclic frequency for an extended 
fluid object Kr as 



2e-^^n [2n + 'uj'^^ 



— 2zun^ 

dzu dzu 



zu dP 
E + Pd^ 



(43) 

A number of comments should now be made about this ex- 
pression. Firstly, while the definition <43> reduces to the corre- 
sponding expression <20> when the Newtonian limit is taken, it has 
a new important correction coming from pressure gradients (cf . sec- 
ond term in the second round brackets, which is a function of zu 
only). Secondly, in those cases in which the (radial) pressure gra- 
dients can be neglected (e.g. for very thin discs), the fluid motion 
is essentially Keplerian and expression 143 > coincides with the rel- 
ativistic radial epicyclic frequency for a point-like particle. In this 
case, and only in this case, <43> reduces to the expression for the 
radial epicyclic frequency derived by Okazaki et al. (1987) that, for 
a Schwarzschild spacetime in spherical coordinates, is 



(Kr) 



Keplerian 



6M \ GM 



(44) 
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Model 


^-distribution 


L 


a, (3 


q 




^out 




"^max 






(cl) 


const. 


3.90 






4.241 


36.26 


4.233 


9.457 


-1.0 X 10- 


-6 


(c2) 


const. 


3.75 






4.882 


11.57 


4.836 


7.720 


-1.0 X 10" 


-5 


(c3) 


const. 


3.70 






5.538 


8.227 


5.266 


6.921 


-1.0 X 10- 


-4 


(11) 


linear 




0.03, 3.70 




5.115 


46.71 


4.444 


10.13 


-1.0 X 10- 


-4 


(12) 


linear 




0.01, 3.70 




5.224 


12.74 


4.845 


8.087 


-1.0 X 10- 


-5 


(pll) 


power-law 


3.70 




0.02 


4.872 


20.97 


4.520 


9.108 


-1.0 X 10- 


-5 


(pl2) 


power-law 


3.60 




0.05 


5.065 


61.98 


4.312 


11.06 


-1.0 X 10- 


-5 


(pl3) 


power-law 


3.66 




0.02 


5.157 


15.46 


4.694 


8.596 


-1.0 X 10- 


-5 


(pl4) 


power-law 


3.68 




0.03 


4.854 


32.13 


4.395 


9.889 


-1.0 X 10- 


-6 



Table 1. Main properties of the two non constant angular momentum models considered. From left to right 
the columns report: the type of angular momentum distribution (see SectionlTlfor a description of £c,a, 13 
and q), the inner and the outer radii of the torus roin and roout, the radial position of the maximum density 
in the torus romax, and the value of the potential gap chosen at the inner edge of the torus. 



An important feature of the relativistic radial epicyclic fre- 
quency j44t . and that distinguishes it from the corresponding ex- 
pression <43t . is that such a function is not monotone in zu but 
has a maximum at a given radius, thus indicating that in General 
Relativity oscillation modes can be trapped near the inner edge of 
an accretion disc (Kato and Fukue, 1980). When pressure gradi- 
ents are taken into account, however, this ceases to be true and the 
epicyclic frequency becomes monotonically decreasing with radius 
(cf. Figure|5j 

Finally, it should be noted that, as for the Newtonian case, 
expression j43> predicts a zero epicyclic frequency for orbital mo- 
tion with constant specific angular momentum. Note that this holds 
true only when the specific angular momentum is defined as £ = 
—u^/ut, in which case 



M d , 2 
avj avj 



'Q,) — zoe 



20. + TO — 

dvj 



dv 

dzu 



(45) 

To the best of our knowledge, this result has not been dis- 
cussed before in the literature for finite-size fluid configurations. 
Furthermore, it is relevant for the physical interpretation of the 
modes found in the numerical calculations of Zanotti et al. (2003) 
as it indicates that the oscillations found in those simulations cannot 
be associated (at least within the validity of a vertically integrated 
approach) with epicyclic oscillations, but should be associated with 
other restoring forces, most notably pressure gradients. 



5 PERTURBATIONS OF RELATIVISTIC TORI: A 
GLOBAL ANALYSIS 

A global analysis of the axisymmetric oscillations modes of a rel- 
ativistic torus consists of solving the system of equations <35t - 
<37t as an eigenvalue problem, treating the perturbed quantities as 
eigenfunctions and the unknown a as the eigenfrequency (see Ro- 
driguez et al., 2002 for the solution of the equivalent problem for 
relativistic thin discs). 

It is worth remarking that equation <28> basically states that a 
condition for the existence of stationary, geometrically thick fluid 
configurations orbiting around a Schwarzschild black hole is that 
the distribution of angular momentum is a non-Keplerian one. As 
long as this condition is met, in fact, the left-hand-side of equation 
<28> will not be zero and the pressure gradients will expand the 
disc in the vertical direction. Of course, different distributions of 



angular momentum will define tori with different properties such 
as: different inner and outer radii, different positions for the maxi- 
mum in density, etc. Among the infinite number of possible angular 
momentum distributions, those that produce tori with a finite radial 
extension are of particular interest, since they allow for the exis- 
tence of globally trapped, axisymmetric modes of oscillation. 

Because of the degeneracy in the functional form for £{taj), 
we have constructed several different models that we have reported 
in Table 1 and that we will be discussing in detail in the following 
Sections. This approach to the problem is clearly more expensive 
as the eigenvalue problem needs to be solved for a large number of 
models; however it also allows one to investigate how the axisym- 
metric oscillations depend on the angular momentum distribution. 
When features are found which do not depend on a specific choice 
for the distribution of angular momentum and that are thus indica- 
tion of a universal behaviour, this approach can be very rewarding 
as we will discuss in the following Sections. 

However, before discussing the results of the numerical so- 
lution of the eigenvalue problems, it is useful to review briefly the 
numerical approach followed in the solution of the eigenvalue prob- 
lem and the boundary conditions imposed. 



5.1 Numerical method 

Equations <35> - <37> . supplemented with suitable boundary condi- 
tions at both edges of the numerical grid, represent a standard two- 
point boundary value problem, for which a "shooting" method can 
be used. In practice, given a trial value for the eigenfrequency a, 
two solutions for the unknown quantity (either SQ or 5U) are found 
starting from the two edges of the numerical grid at the inner, zuin, 
and outer, zuout, radii of the torus. These two solutions are then 
matched at an arbitrary point 7i7„ in the domain, where the Wron- 
skian of the left and right solutions is evaluated. This procedure 
is iterated until a zero of the Wronskian is found with the desired 
accuracy, thus determining both the eigenfrequency and the eigen- 
functions of the mode considered. 

Note that because of the linearity of the present analysis, the 
solutions are invariant after a rescaling of all the quantities by a 
constant. In practice, this freedom allows to choose the initial guess 
for one of the eigenfunctions at one edge of the numerical grid 
to be a freely specifiable number, which we typically set to one. 
Furthermore, since the spacetime is uniquely described once the 
mass Al of the black hole has been fixed, we can use it to rescale 
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all of the relevant quantites involved in the calculations. The new 
dimensionless quantities, indicated with a "tilde", are then defined 
as 



= ZJ 



M \ (GMp 



MqJ \ c2 



E = E 



(^0) ( 



Ma 



P = P 



M_\'^ ( GM, 



® \ Ayr „2 

Mq c , 



(46) 
(47) 
(48) 



K = K 



GMp 



{.Me J y c2 



2{r-i) 



i , (GMe 



{Me)'-^ c , (49) 



(50) 



MqJ\ 

Unless specified differently, all of the results presented hereafter 
will be given in terms of these dimensionless quantities. 

5.2 Boundary conditions 

Suitable boundary conditions need to be specified at the inner and 
outer edges of the vertically integrated torus for the solution of the 
set of equations <35t - <37> . As for oscillations in stars, we here as- 
sume the perturbed surface of the torus to be where the Lagrangian 
pressure perturbation vanishes, i.e. where 

Ap = . (52) 

Because the Eulerian and Lagrangian descriptions are linked 
through the relation 

A = S + C^, (53) 

where is the Lie derivative along the Lagrangian displacement 
three-vector ^, condition <52> can also be written as 

5p + ^'djp^0, (54) 

from which it follows that Ap = — Sp for a polytropic fluid 
configuration whose rest-mass density vanishes at the surface (Tas- 
soul, 1978). Note also that while in the Cowling approximation the 
Eulerian perturbation of the metric, Sgab is zero, the Lagrangian 
perturbation Agab = Sgab + £-^ab = Va^fc + Va^fc is, in general, 
nonzero. 

Using the relation between the perturbed 4-velocity and 
the Lagrangian perturbation of the metric (see Friedman 1978), 
1 
2 

which, in our case, reduces to the relation, 



(51) 



Au" = ^u°'u'^u'Agi3j , 



(55) 



r = 

u 

we can write the boundary condition j52t simply as 



dP 



517 = 0, 



(56) 



(57) 



cr{E + P) duj 

Expression <57L taken as a linear relation between SQ and SU, 
provides us with the boundary conditions to be imposed at the inner 
and outer edges of the torus. 



6 CONSTANT SPECIFIC ANGULAR MOMENTUM TORI 

The simplest and most studied choice for the distribution of angu- 
lar momentum is the one in which i = £c — const., with £c being 
chosen in the interval between the specific angular momentum at 
the marginally stable orbit Ims = 3\/6/2 — 3.67 and the one at 
the marginally bound orbit £mh = 4, to yield a disc of finite size. 
Besides yielding a simpler approach, constant specific angular mo- 
mentum tori benefit from having a radial epicyclic frequency which 
is identically zero [cf. eqs. <43> and <45H . thus leaving pressure gra- 
dients as the only possible restoring forces. As a result, a constant 
angular momentum distribution provides an interesting limit for the 
properties of p-mode oscillations that will not be influenced by cen- 
trifugal effects. 

The procedure for the construction of the background model 
for the torus can be summarized as follows. Firstly, after a value 
£c for the constant specific angular momentum has been chosen , 
the locations of the cusp rocusp and of the maximum density in the 
torus^ ti7niax are determined as the positions at which the Keple- 
rian specific angular momentum is equal to £c. The inner and outer 
edges of the torus can then be calculated after requiring that they 
lie on the same equipotential surface, i.e. 



Ut{uJin) = Util^cuBp) exp(AVKin) = Ut{uJout) 



(58) 



where AWin is the value of the potential gap chosen at the inner 
edge of the torus and is therefore a measure of how much the torus 
"fills" its Roche lobe. Clearly, the case in which AWin ~ corre- 
sponds to a torus filling its outermost closed equipotential surface 
(i.e. the one possessing a cusp) 

Having determined vjin, it is then possible to calculate the an- 
gular velocity distribution as 



2M \ 
- 2M ) 



(59) 



where Qc is angular velocity at the radial position of cusp cc7cusp. 
Finally, after making a choice for the parameters of the equation 
of state (and therefore for k and 7) the equation for the hydrostatic 
equilibrium <31> can be integrated, either numerically or analyti- 
cally (Kozlowski et al., 1978), to build the torus model which will 
serve as background for the introduction of the perturbations. 

When the specific angular momentum is constant within the 
torus, also the perturbation equations become simpler and, in par- 
ticular, the second term in <36> vanishes exactly [cf. equation <45H 
and the resulting set of equations can then be written as 



aSU + Ae 



. djSQ) 
dvj 



■ 



5V^ + AQ.e''^"vjSQ = , 



(60) 



(61) 



a5Q + e" 



VP djSU) 

^ du 
dzu 



P dvj 
1 



E + P 



1 dA ^ 1 dP 
A dzu rP dzu 



3U = 0. (62) 



Equations <60> - <62t have been solved numerically for a num- 
ber of different models, whose main properties, such as the inner 



^ Note that the position of the maximum density in the torus is traditionally 
refen'ed to as the "centre" of the torus, although its coordinate centre is, of 



course, at ro : 
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Figure 1. Eigenfunctions for 5Q = 5P/ {E + P) as a function of the 
radial coordinate for a constant angulai' momentum disc. Only the funda- 
mental mode / and the first three overtones, denoted by oi, 02 and 03, 
have been reported. The data refers to model (c2) of Table 1 and the units 
on the vertical axis are arbitrary. 



Figure 2. Eigenfunctions for 5U = iSV^ as a function of the radial coor- 
dinate for a constant angular momentum disc. Only the fundamental mode 
/ and the first three overtones, denoted by oi, 02 and 03, have been re- 
ported. The data refers to model (c2) of Table 1 and the units on the vertical 
axis are arbitrary. 



and outer radii roin, Wout, tlie position of the maximum in the sur- 
face density distribution Sjmax, and the value of the specific an- 
gular momentum at the inner edge of the disc £c = £{u:'in), are 
summarized in Table 1 . As a representative example, we present in 
FiguresQand 2 the results for a torus model with £c = 3.75 [i.e. 
model (c2) of Table 1]. 

Figure Q in particular, shows the solution for the eigenfunc- 
tion SQ for both the fundamental mode / and for the first three 
overtones 01,02 and 03. Similarly, Figure|2|shows the fundamental 
mode and the first three overtones for the eigenfunction SU, which 
has always one node less than 5Q. Note that the choice of a very 
small potential barrier introduces complications in the numerical 
solution of the eigenvalue problem. For a model with AWin — 0, 
in fact, the fluid elements at the inner edge are just marginally sta- 
ble to accretion onto the black hole. This behaviour is reflected in 
the eigenfunctions SQ and 5U, which tend to diverge at the inner 
and outer edges of the torus. To avoid this problem (still partially 
visible in FiguresQand|2} we have used a small but nonzero value, 
AWin = — 10~^, for the potential barrier. 

A careful look at FiguresQand|2|reveals that the modes / and 
02 have similar properties at the inner edge; furthermore, during 
each half-period, they have signs which are opposite to those of the 
modes oi and 03. This is an important feature indicating that not 
all modes will behave in the same manner at the inner edge where 
mass-loss and accretion onto the black hole take place. While a con- 
clusion on the this behaviour cannot be drawn on the basis of the 
present linear analysis, in which the total mass is conserved exactly 
(E5U = at the inner and outer edges for all modes), the func- 
tional form of the eigenfunctions for 5Q and SU seems to suggest 
that mass accretion through the cusp will be possible preferably at 
the frequencies corresponding to the /, and 02n modes (see also the 



discussion below for a comparison with the numerical calculations 
presented in Zanotti et al., 2003). 

The eigenfrequencies corresponding to model (c2), as well as 
for all of the other models in Table 1, are listed in Table 2, where 
we have included the frequencies up to the fourth overtone. A care- 
ful look at Table 2 reveals that the computed eigenfrequencies are, 
at least for the first few modes, in a sequence 2:3:4:. . ., to a good 
precision. As the order of the mode increases, this simple sequence 
is no longer followed and a different one appears. The fact that the 
lowest order modes appear at frequencies that are in small integer 
ratios is not particularly surprising if p modes are to be interpreted 
as sound waves trapped within the cavity represented by the con- 
fined torus. In this case, in fact, one would indeed expect that modes 
should be trapped and with wavelengths that are multiples of the 
trapping lengthscale, i.e. A„ — [2/(2 + n)]L, where n = 0, 1, . . .. 
Stated differently, the results reported in Table 2 are consistent with 
the idea that p modes can be associated with sound waves trapped 
in the torus and having frequencies which are multiples of the fun- 
damental one and given by o-„ = Cs/X„ = [(2 + n)/2]f. 

Note that while these p modes follow the same definition and 
are similar in nature to the ones discussed by Nowak and Wag- 
oner (1991, 1992) for thin discs, their trapping is not produced by 
a turning of the relativistic radial epicyclic frequency at smaller 
radii. Rather, the relativistic radial epicyclic frequency in discs in 
which pressure gradients play an important role, has been found to 
be monotonically decreasing for all of the disc models considered 
here [cf. equation <43l , Fig.|5)- An important consequence of this 
is that the p modes for thick discs are not restricted to be trapped 
in special parts of the disc, as it happens for the corresponding p 
modes for thin discs (Nowak & Wagoner 1991, 1992). In particu- 
lar, they need to be present only in the (rather small) inner regions 
of the disc, where they can produce only a slight modulation in the 
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Model 


/ 


Ol 


02 


03 


04 


(cl) 


0.00838 


0.01227 


0.01604 


0.01973 


0.02337 


(c2) 


0.02017 


0.02950 


0.03700 


0.04322 


0.04869 


(c3) 


0.01898 


0.02764 


0.03514 


0.04230 


0.04933 


(11) 


0.00602 


0.00831 


0.01060 


0.01287 


0.01514 


(12) 


0.01992 


0.02856 


0.03618 


0.04330 


0.05019 


(pU) 


0.01447 


0.02114 


0.02733 


0.03324 


0.03894 


(pl2) 


0.00426 


0.00616 


0.00803 


0.00988 


0.01172 


(pl3) 


0.01810 


0.02635 


0.03764 


0.04073 


0.04742 


(pl4) 


0.00948 


0.01379 


0.01792 


0.02196 


0.02593 



Table 2. Eigenfrequencies of the fundamental / mode and of the first four 
overtones o„ of the vertically integrated models described in Table 1 . All 
the frequencies are given in normalized units. 



disc emission. Similarly, they are not constrained to be trapped the 
outer regions of the disc, where the location of the outer radius re- 
mains uncertain and the corresponding frequencies very small. On 
the contrary, the modes discussed here are present over the whole 
torus, which then behaves as a single trapping cavity. This is a sub- 
stantial difference, giving these modes the possibility of possess- 
ing frequencies comparable with observations and of causing non- 
negligible modulations in the disc emission. 

Having solved the eigenvalue problem for p modes in ver- 
tically integrated relativistic tori with constant specific angular 
momentum, we can now proceed to a comparison with the fully 
nonlinear two-dimensional simulations discussed by Zanotti et al. 
(2003). This is presented in Figure |3| where we have plotted the 
power spectra of the L2 norms of the rest-mass density for the torus 
model (c2) in Table 1. The solid line, in particular, refers to initial 
data in which a global perturbation has been introduced in terms of 
the radial velocity field of a relativistic spherical accretion solution 
[cf. equation (15) of Zanotti et al., 2003]. Clearly, the power spec- 
trum for the continuous line shows the presence of peaks appearing 
in a simple small integer sequence 2:3:4:. . ., providing convinc- 
ing evidence that the oscillations triggered in the simulations of 
Zanotti et al. (2003) correspond indeed to p-modes oscillations of 
perturbed relativistic tori. It should also be noted that this sequence 
is not the same shown by the peaks in the power spectrum of the 
mass accretion rate onto the black hole, which instead are only in 
the sequence 1:2:3:. . . (cf. Figure 7 of Zanotti et al., 2003). This 
difference remains puzzling, but could be explained on the basis 
of the behaviour of the eigenfunctions for the / and 02n modes 
at the inner edge discussed above. In this case, in fact, all of the 
02n+i modes would have vanishingly small perturbations at the in- 
ner edge of the disc and would not produce a mass accretion at 
those frequencies, as shown by Figure 7 of Zanotti et al. (2003). 

The knowledge of the eigenfunctions of p modes in relativis- 
tic tori can also be used to go a step beyond a simple comparison 
with numerical simulations and actually use the latter to perform in- 
vestigations of "numerical" discoseismology. In other words, hav- 
ing a complete picture of the oscillation properties and an accurate 
two-dimensional numerical code to simulate them, it is possible to 
investigate the dynamical response of a relativistic thick disc as a 
result of the introduction of suitably selected perturbations. As a 
concrete example, we report with the dashed line in Figure |3| the 
power spectrum of the L2 norm of the rest-mass density for the 
torus model (c2) in Table 1 . The important difference with the cor- 
responding spectrum indicated with a solid line is that the dashed 
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Figure 3. Power spectra of the L2 norm of the rest-mass density for the 
torus model (c2) in the Table 1 of Zanotti et al. (2003). Different lines refer 
to models with different initial perturbations. The soUd line, in particular, 
refers to initial data in which a global perturbation has been introduced in 
the radial velocity field only [cf. equation (15) of Zanotti et al. (2003)]. The 
dotted line, on the other hand, refers to initial data in which a perturbation 
based on the computed eigenfunctions of the oi mode has been introduced 
both in the radial velocity and in the density. The two spectra have been 
reseated in order to match in the power of the fundamental frequency. 



line refers to a simulation having as initial data perturbations based 
on the computed eigenfunctions for 5U and SQ of the oi mode. 
This represents a selective excitation of the oi mode and, as a re- 
sult, the corresponding power in the first overtone is increased by 
almost a factor of ten (the power in the 02 mode is instead de- 
creased of a similar amount). This behaviour further confirms that 
the p modes derived in this perturbative analysis correspond to the 
modes simulated numerically by Zanotti et al. (2003). Finally, it 
should be noted that the power spectrum corresponding to initial 
data containing only an oi-mode perturbation shows that also the 
other modes have been excited, most notably the fundamental one. 
This is due partly to a mode-mode coupling which transfers energy 
from one mode to the other ones, but also to the error introduced 
using eigenfunctions derived for a vertically integrated model, and 
that cannot reproduce the corresponding eigenfunctions of a two- 
dimensional disc exactly. 

There is a final aspect of the axisymmetric p modes dis- 
cussed so far which is worth underlining. This is illustrated in Fig- 
ure|4]where we have plotted the values of the eigenfrequencies for 
the fundamental mode / (solid line) and for the first overtone oi 
(dashed line) for large number of tori. All of the points on the solid 
and dashed curves in Figure |4] represent the solution of the eigen- 
value problem and the two sequences have been calculated for tori 
with fixed "centre" t&max = 10.17, but having different radial ex- 
tensions L = 57out — fijin. As one would expect for modes behaving 
effectively as sound waves, the eigenfrequencies decrease as the ra- 
dial extent of the torus increases. Furthermore, the frequencies for 
/ and Ol maintain a harmonic ratio 2:3 for all the values of L. Most 
importantly, however, the filled dot shown in Figure |4| for L — 
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Figure 4. Eigenfrequencies for the fundamental / (solid line) and for the 
first overtone o\ (dashed line) of axisymmetric p modes for i = const, 
tori. The two lines refer to sequences of tori having the same radial po- 
sition for the "centre" romax = 10.17. but different radial extensions 
L = roout — roin- The filled dot represents the radial epicyclic frequency 
for a circular orbit at romax and thus the value at which the fundamental 
p-mode frequency tends in the limit of a vanishing torus size. 



represents the radial epicyclic frequency for a circular orbit at Sjc, 
i.e. a = ^(1 - &/^c)GM/vj'i. 

This is an important result that provides a simple interpreta- 
tion of the nature of the axisymmetric p modes of thick discs: the 
radial epicyclic frequency at the position of maximum density rep- 
resents the value at which the fundamental p-mode frequency tends 
in the limit of a vanishing torus size. The importance of this result 
is that while different models for the tori, using, for instance, differ- 
ent equations of state or different angular momentum distributions, 
will produce different slopes for the solid and dashed lines in Fig- 
ure |4] all of the sequences will terminate on the filled dot when 
L — > 0, that is, when finite-size effects will cease to be relevant 
and each torus will effectively behave as a ring of particle in circu- 
lar orbits. Because this occurs only in the limit L ^ 0, it should 
not surprise that it is valid also for sequences of tori with constant 
distributions of specific angular momentum for which, as discussed 
in Section|4] the epicyclic frequency is zero. Indeed, this result is 
totally independent of the distribution of specific angular momen- 
tum and will be confirmed in the following Section, where tori with 
non-constant specific angular momentum distributions are consid- 
ered. An analytic proof of this conclusion is provided in Appendix 
A. 



7 NON-CONSTANT SPECIFIC ANGULAR MOMENTUM 
TORI 

While the study of constant angular momentum discs offers several 
advantages and simplifies the equations, realistic disc are likely to 
have angular momentum distributions that are not constant. For this 
reason, and in order to assess the validity of the results derived so 



far with more generic distributions of angular momentum, we have 
extended our mode analysis also to discs with non-constant angu- 
lar momentum. The first step in this direction consists of determin- 
ing which distribution of specific angular momentum I = i[-uj) 
should be specified in the construction of an equilibrium model. 
This choice is, to some extent, arbitrary with the only constraint be- 
ing given by Rayleigh 's criterion for the dynamical stability against 
axisymmetric perturbations (Tassoul 1978). This condition, how- 
ever, is not very strong and simply requires that di/dzo > 0, so 
that even a constant angular momentum distribution is stable, al- 
though only marginally. 

Some guidance in this choice can come from numerical sim- 
ulations and indeed calculations of torus formation performed by 
Davies et al. (1994) with smooth particle hydrodynamics (SPH) 
techniques have shown that the final configuration consists of a core 
object surrounded by a torus whose angular momentum distribution 
is close to a power-law in radius, with index 0.2. In our computa- 
tions we have therefore adopted two simple and different distri- 
butions for the specific angular momentum that can be expressed 
analytically. We refer to the first one as the "linear" distribution, in 
which 



OIZU + f3 . 



(63) 



where a and /3 are adjustable constant coefficients. Similarly, we 
refer to the second one as the "power-law" distribution, in which 
instead 



(64) 



where q is also an arbitrary parameter. The values chosen for 
a, 13, Ic and q for the representative models discussed here have 
been summarized in Table 1 . Note that we have deliberately chosen 
small values for both a and q and this is to avoid the construction of 
unpertubed models that differ significantly in radial extension from 
the ones built with a constant specific angular momentum distribu- 
tion^. 

The construction of the equilibrium models proceeds in this 
case as discussed in the previous Section, namely, by first fixing an 
initial value for the specific angular momentum at the cusp, then by 
calculating the radial extension of the torus in terms of the potential 
barrier AWin, and finally by integrating numerically equation M9\ 
for the chosen distribution of Q.{-uj). 

The full system of the perturbed equations <35t - <37t . can now 
be rewritten in a more compact form as 
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2e-nfl + ^£).H., (66) 



= a5W+{ ^ + -n-2n^]u,e-''SU , (67) 
dzu zu dzu , 



It is not difficult to realize that a constant specific angular momentum dis- 
tribution yields the most compact tori among those having the same angular 



momentum at the cusp. 
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Figure 5. The upper panel shows the relativistic sound velocity for rep- 
resentative tori models having either a constant [solid line; model (c2)], 
a linear [short-dashed line; model (12)] or a power-law [long-dashed line; 
model (pl4)] distribution of specific angular momentum. The lower panel, 
instead, shows the dimensionless radial epicychc frequency for the three 
models shown in the upper panel, with the horizontal lines showing the fun- 
damental frequencies for models (12) and (pl4). All quantities are expressed 
in normalized units. 



and can be solved adopting the same numerical technique used for 
the constant angular momentum models. 

The upper panel of Figure|5|shows the relativistic sound veloc- 
ity for representative unperturbed tori models having either a con- 
stant (solid line), a linear (dotted line) or a power-law (dashed line) 
distribution of specific angular momentum [The data refers to mod- 
els (c2), (12) and (pl4) of Table 1, respectively.]. Note that while all 
of the curves refer to models with a polytropic EOS, each curve 
does not depend on the value chosen for the polytropic constant k 
and hence on the mass of the disc (see Appendix B for a proof of 
this). The lower panel of Figure|5| on the other hand, shows the ra- 
dial epicyclic frequency Kr as calculated from expression <43j for 
the three models shown in the upper panel (we recall that = 
when £ — const.). The horizontal lines in the lower panel refer 
to the fundamental frequencies for models (12) and (pl4) and their 
relative position with respect to the curves for Kr will be important 
for the appearence of an evanescent-wave region in the inner parts 
of the torus (see the discussion below). 

We have summarized in the different panels of Figure Q the 
eigenfunctions for 5Q, SU and SW for tori with a linear [panels 
(a)-(c)] and a power-law [panels (d)-(f)] distribution of specific an- 
gular momentum. The eigenfunctions refer, in particular, to models 
(11) and (pl3) of Table 1 and are shown in the fundamental as well 
in the first three overtones. Similarly, are summarized in Table 2 
also the eigenfrequencies computed for a number of tori with non- 
constant specific angular momentum. 

A rapid analysis of the eigenfrequencies in Table 2 as well as 
of the eigenfunctions in the various panels of Figure Q indicates 
that, overall, the results for non-constant distributions of specific 
angular momentum do not differ, at least qualitatively, from those 
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Figure 6. Schematic propagation diagram for a thick disc. The thick soHd 
hne shows the values of the relativistic radial epicyclic frequency Kr in the 
inner regions of the disc, while the horizontal dashed line the value of the 
fundamental frequency /. The data is presented in normalized units and 
refers to model (pl2) of Table 1 . See the main text for a discussion. 



computed for constant distributions. Most importantly, the eigen- 
frequencies continue to appear in the simple sequence 2:3:4:. . ., at 
least for the lowest-order modes considered here. 

This result may appear surprising in light of the fact that the 
the radial epicyclic frequency defined by <43t does not vanish in 
the case of non-costant specific angular momentum distributions, 
so that the computed eigenfrequencies contain also a contribution 
coming from the inertial oscillations [cf. the dispersion relation 
J42H . However, it is important to bear in mind that the p modes 
discussed here behave essentially as sound waves trapped inside 
the torus. As a result, while their absolute frequencies will be dif- 
ferent when tori with different radial extensions are considered, the 
sequence in which they appear will basically follow the relation 
o-„ = [(2 + n)/2]/. 

The "universality" in which the lowest-order p-mode seem 
to appear has an interest of its own, but may also have important 
implications for the high frequency quasi-periodic oscillations ob- 
served in binaries containing a black hole candidate (Strohmayer 
2001; Remillard et al., 2002). In at least three systems, in fact, the 
QPOs have been detected with relatively strong peaks that are in 
a harmonic ratio of small integers 1:2, 2:3, or 1:2:3 (Abramowicz, 
Kluzniak 2001), and a possible explanation of this phenomenology 
in terms of the oscillations of a geometrically thick, non-Keplerian 
disc of small radial extent orbits around the black hole will be pre- 
sented in a separate paper (Rezzolla et al., 2003). 

Before concluding this Section, we should underline an im- 
portant qualitative difference that emerges when non-constant dis- 
tributions of specific angular momentum are considered. This is 
related with the behaviour of the eigenfunctions for the radial ve- 
locity perturbation SU at the inner edge of the disc. A rapid com- 
parison of Figure |2| with the corresponding panels (b) and (e) of 
Figure shows that the eigenfunctions become vanishingly small 
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Figure 7. Eigenfunctions for SQ = 5P/{E + P), SU = iSV^, and SW = SV^ as a function of the radial coordinate for a tori with a "hnear" [panels 
(aMc); model (11)] and a "power-law" [panels (d)-(f); model (pl3)] distribution of the specific angular momentum. While the units used on the vertical axes 
are arbitrary, it should be noted that the values on panel (e) are one order of magnitude smaller than the corresponding ones on panel (b). In all panels only the 
fundamental mode / and the first three overtones have been reported. 
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at the inner edge of the disc in the case of "power-law" distribution 
of specific angular momentum. This does not seem to be case for 
the behaviour at the outer edge, where the eigenfunctions maintain 
very large values. Similarly, this behaviour at the inner edge is not 
encountered with the other two distributions considered [note that 
the values on vertical axis of panel (e) are one order of magnitude 
smaller than the corresponding ones on panel (b) and in Figure|2|. 

A physical explanation for this behaviour can be found by 
comparing the behaviour of the radial epicyclic frequency curves 
for the different distributions of angular momentum. The lower 
panel of Figure |3| in particular, shows that the radial epicyclic fre- 
quency for both the linear and the "power-law" distributions of 
specific angular momentum are monotonically decreasing for in- 
creasing radii (This has been found to hold for all of the disc mod- 
els considered here.). However, while the two frequencies tend to 
be comparable in the outer regions of the torus, they differ con- 
siderably in the inner regions, where the epicyclic frequency for a 
power-law distribution increases inward much more rapidly. More 
importantly, the epicyclic frequency is always smaller than the fun- 
damental frequency for a torus with a linear distribution of angular 
momentum, while this is not the case for a power-law distribution 
(cf. Table|2j. 

In this latter case, then, there will be regions of the disc in 
which cr^ — ~ fc^Cs < and thus where the amplitude of 
any acoustic wave is forced to be zero (Other modes, such as the g 
modes, can however exist in this region, in which they are actually 
trapped; Kato and Fukue, 1980). Such a region is referred to as an 
evanescent-wave region and effectively represents the part of the 
disc in which the the centrifugal barrier produced by rotation makes 
the radial epicyclic frequency so high that all inwarding moving 
acoustic waves cannot propagate but are reflected back to larger 
radial positions. In those parts of the disc where the fundamental 
frequency is larger than the radial epicyclic one, on the other hand, 
a^ — K^ > 0, and sound waves can propagate undamped at least in a 
perfect fluid. Stated differently, the results found here indicate that, 
when acoustic waves are considered, the inner regions of the torus 
are of evanescent type for discs with a power-law distribution of 
angular momentum, while no evanescent region exists in a constant 
or in a linear distribution of specific angular momentum. 

Finally, it is worth pointing out that the behaviour discussed 
above can be used to interpret some of the most recent results on 
the occurrence of the runaway instability (Font and Daigne 2002a, 
2002b; Zanotti et al., 2003). In a series of relativistic hydrodynam- 
ics calculations of a torus orbiting and accreting onto a rotating 
black hole, in fact. Font and Daigne (2002b) have found that while 
tori with constant specific angular momentum distributions lead to 
a runaway instability, a slight outward increase as a power-law of 
the specific angular momentum, has a dramatic stabilizing effect, 
suppressing the instability completely. While this behaviour reflects 
a process which is fundamentally nonlinear, it is indeed consistent 
with the interpretation that the inner regions of thick discs prevent 
the inward propagation of perturbations and hence tend to suppress 
the accretion of mass onto the black hole. Because the runaway 
instability is fed by the changes in the spacetime produced by the 
mass accretion, a considerable reduction of the latter can be the 
cause of the suppression of the former. A more detailed investiga- 
tion is necessary to validate this hypothesis. 



8 CONCLUSIONS 

We have presented the first investigations of the oscillation proper- 
ties of relativistic, non- self gravitating tori orbiting around a black 
hole. More precisely, we have here considered the axisymmetric 
oscillation modes of a geometrically thick disc constructed in a 
Schwarzschild spacetime. Like relativistic stars, relativistic tori are 
extended objects with non-trivial equilibrium configurations gov- 
erned by the balance among gravitational, centrifugal and pressure 
forces. Unlike relativistic stars, however, the contributions coming 
from centrifugal forces are not small in relativistic tori and this am- 
plifies the role of oscillations restored by centrifugal forces. This 
leads to oscillations properties that can be considerably different 
from those encountered in stars. 

Because thick discs have angular momentum distributions that 
are intrisically non-Keplerian, we have modeled them with a num- 
ber of different distributions of specific angular momentum, which 
have been chosen to be either constant within the torus, or with a 
radial dependenced being linear or a power-law. Furthermore, in 
order to keep the treatment as simple as possible and handle equa- 
tions analytically whenever possible, we have built the models for 
the tori using vertically integrated and vertically averaged quanti- 
ties. 

Rather little is still known about the oscillation modes of rel- 
ativistic tori and our approach to the problem has therefore pro- 
gressed by steps. In particular, our first step has been that of consid- 
ering the dispersion relation for acoustic waves propagating within 
these objects. We have done this firstly in a Newtowian framework 
in which equations are simpler and so is their physical interpreta- 
tion. We have then extended the local analysis to a general rela- 
tivistic framework and in particular to a Schwarzschild black hole 
spacetime. While a local analysis and the resulting dispersion re- 
lation is valid only for oscillations with wavelengths small when 
compared with the typical lengthscale in the tori, the ones derived 
here have been important to clarify the relation between the acous- 
tic waves and the other waves that play a fundamental role in fluids 
orbiting in a central potential, i.e. inertial (or epicyclic) waves. In 
particular, it has been possible to show that, in general, both acous- 
tic and inertial oscillations are present in the perturbations of thick 
discs and that the latter can be removed only in the special case 
of a constant distribution of specific angular momentum. This re- 
sult, which was already know in Newtonian physics, has here been 
shown to hold also for extended relativistic fluid configurations. To 
to best of our knowledge this result has never been discussed in the 
literature before. 

Going a step further and beyond a local analysis, we have used 
the same mathematical setup to perform a global analysis and de- 
termine both the eigenfunctions and the eigenfrequencies of the ax- 
isymmetric oscillations. Also in this case, the assumption of a verti- 
cally integrated equilibrium has simplified the eigenvalue problem 
considerably, translating it into a set of coupled ordinary differen- 
tial equations which have been solved using standard techniques 
and for a number of different models for the tori. The modes found 
in this way correspond to the p modes of relativistic tori and are 
characterized by eigenfrequencies appearing in a simple sequence 
of small integers 2:3:4:. . ., at least for the first lower-order modes. 
This important feature does not depend on the distribution of angu- 
lar momentum used to build the tori and basically reflects the fact 
that p modes are to be interpreted as sound waves trapped within 
the cavity represented by the confined torus. 

The properties found here with a linear global analysis are in 
very good agreement with the numerical results found in the time 
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evolution of perturbed "toroidal neutron stars" having constant dis- 
tributions of specific angular momentum (Zanotti et al., 2003). This 
agreement with fully nonlinear two-dimensional simulations pro- 
vides convincing evidence that the assumption of a vertically inte- 
grated equilibrium has not subtracted important information from 
the present analysis. 

A number of motivations are behind the investigations carried 
in this paper. Firstly, and as mentioned above, relativistic tori pos- 
sess a rich variety of oscillation modes which are still essentially 
unexplored, in stark contrast with what is known (both in New- 
tonian and in relativistic regimes) about geometrically thin discs. 
In this respect, the present work was intended as a first investiga- 
tion of the discoseismology of geometrically thick discs. Secondly, 
our analysis was aimed at interpreting and clarifying some of the 
results found by Zanotti et al. (2003), who had investigated numer- 
ically the dynamical response of massive tori to perturbations but 
had not been able to characterize all of the mode properties. Finally, 
by investigating the response of tori to perturbations, we intended 
to highlight possible connections between the oscillation modes of 
these objects and all those astrophysical scenarios in which geo- 
metrically thick discs may play an important role. 

Some of the results discussed here could find application in 
more general contexts. One of such applications is offered by 
the quasi-periodic X-ray phenomenology observed in X-ray bi- 
nary systems containing a black hole candidate (Strohmayer 2001; 
Remillard et al., 2002). In these systems, in fact, the X-ray emission 
is modulated quasi-periodically with power spectra having peaks 
in a harmonic ratio of small integers 1:2, 2:3, or 1:2:3 (Abramow- 
icz and Kluzniak, 2001). While there are a number of possible ex- 
plantions for this behaviour, it could be interpreted simply in terms 
of the p-mode oscillations of a geometrically thick, non-Keplerian 
disc of small radial extent orbiting around the black hole (Rezzolla 
et al., 2003). Another of such applications is offered by the runaway 
instability which has recently been considered by a number of au- 
thors (Font and Daigne 2002a, 2002b; Zanotti et al., 2003) through 
fully nonlinear numerical calculations. In particular, the numerical 
evidence that the instability can be suppressed if the torus has a 
distribution of angular momentum which is slightly increasing out- 
wards, can be interpreted simply in terms of the behaviour of the 
eigenfunctions for the radial velocity perturbations, which become 
vanishingly small in the inner regions of the disc. This, in turn, is 
produced by the appearence of a region of evanescent-wave prop- 
agation in the inner regions of the torus that prevents the inward 
propagation of perturbations, reduces the accretion of mass onto 
the black hole, thus suppressing the instability. 

As mentioned in the abstract, the present work represents the 
first of a series of papers aimed at a systematic investigation of the 
oscillation properties of relativistic tori. We are presently investi- 
gating extensions of the present approach in which the p modes are 
investigated both when the background spacetime is that of a Kerr 
black hole (Rezzolla & Yoshida, in preparation), and when devi- 
ations from axisimmetry are present. In addition to this, work is 
also in progress to extend the solution of the eigenvalue problem 
to a fully two-dimensional model for the torus. The results of these 
investigations will be presented in future papers. 
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CT-2000-00137). The computations were performed on the Be- 
owulf Cluster for numerical relativity "AlbertlOO", at the Univer- 
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APPENDIX A: ON THE VALUES OF THE 
EIGENFREQUENCIES WHEN L^O 

In the limit of vanishing torus sizes, L ^ 0, all of the eigenfunc- 
tions can be assumed to have only a simple linear dependence on 
zu, i.e. 



5U — a,zu + b . 



SW ^ cuj + d, 



5Q = ezu + f , 



(Al) 



(A2) 



(A3) 



where a, b, . . . , f are just constant coefficients. Using this ansatz in 
the system of perturbed equations <35> - <37> . we obtain the follow- 
ing system of equations 

(jh + 2Ae'^nS{zo) +ikAe~'''^i=0 , (A4) 



■oj\ac + &e ■''//(iaj)] + ad + be '^H[tjj) = Q. 



(A5) 



(re H — ae ik 

E + P 



+ af 



VP 

+ ^ , ^ e^-^ikh = , (A6) 



E + P 



where 



H(vj) = 20. + vj- 2Q.-CO — 

dvj dvj 



S{vj) = 1 + 



VJ dP 
E + Pd^ 



(A7) 
(A8) 



and where it is easy to realize that the radial epicyclic frequency is 
then simply given by [cf. eq. <43H 



= 2e''^^nH{zu)S{zu) 



(A9) 
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Each of the equations in the system <A4t - fA6l can be written 
symbolically as aizu + Pi = 0, i — 1,2,3 with and Pi being 
just shorthands for the longer expressions in <A4t - fASl . Further- 
more, because each equation in of the <A4t - lA6l must hold for any 
uj, we need to impose that all of the coefficients and Pi vanish 
simultaneously. This generates a system of six equations in the six 
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unknowns a, b, . . . , f 

aa. + 2cile~'^ S(taj) + ikAe^'^~^e = 



ah + 2dne~^ S{77^) +ikAe~"'^f = 



ac + ae H{zu) — , 



ad + he~^H{zu) = 0, 
ae + j— ilea — , 

r p 

af+ e'^-^ikh = 0, 

S + P 

which, in matrix form, can be written as 



/ a ^ 

b 



(AlO) 
(All) 
(A12) 
(A 13) 
(A14) 
(A15) 



A 



- . 



(A16) 



V f / 

A non-trivial solution to the system (A16) is possible if the de- 
terminant of the matrix A vanishes and, after lenghty but straight- 
forward calculations, this condition yields 



det{A) = a^ [cr^ - 2fle''^^H{u^)S{z 



V P 

2ne-''HMSiuj)+a(-^-^ 



k^A^e-^' + 



2^^k'Ae~^' 
E + P 



= (A17) 



Since FP/E + P ~ c^, all of the terms proportional to this quan- 
tity can be neglected in <A17t . which then reduces to 

cr^ [cr^ - 2ne''^^H{uj)S{zu)j ^ = . (A18) 
A non-trivial solution of equation <A18> is then given by 



a^ = 2ne-'^^H{uj)S{zu) = k] 



2 

r ? 



(A19) 



thus proving that the eigenfrequencies tend to the value of the 
epicyclic frequency in the limit ot L — * 0. 



The right-hand-side of equation <B1> depends only on the kine- 
matic and geometric properties of the disc, namely on the angu- 
lar momentum distribution (through Q.) and on the external grav- 
itational field of the central black hole (through i/). Both of these 
quantities depend on zu only. 

Introducing the relation 7 = 1 + 1/n, the polytropic EOS 
p — kp^ can also be written in terms of the Emden function 
e" = p as 



N+l 



p = kQ 

so that equation <BH effectively becomes 



1 + V 



where we have defined ijj = k{n+ 1)0. 

Equation jB3> can be integrated analytically to give 



tp{i^) = C exp 



f{zu')dvj' 



(B2) 



(B3) 



(B4) 



where C is an arbitrary constant. 

The important result contained in equation <B4t is that i/) is a 
function of uj only and will not depend, therefore, on the specific 
value chosen for k. As a result, once 7 (and thus n) is fixed, any 
transformation K Ka, where a is an arbitrary constant, must 
be accompanied by a corresponding transformation O Q/a. An 
important consequence of this is that all of the quantities given by 
p/p, dp/dp, p/e, dp/de are invariant under changes in k. Equally 
invariant is the sound speed defined as Cs = \J dp/de. 

To appreciate the physical implications of this result, it is use- 
ful to recall that the polytropic constant plays here the role of deter- 
mining the total rest-mass of the torus. Once the other parameters 
of the background model have been fixed (i.e. the distribution of 
angular momentum and the potential gap), in fact, the rest-mass of 
the torus is calculated from the integral 



2'K'GoTiu''ujdzu 



(B5) 



As a consequence, different choices for k will yield tori with dif- 
ferent density distributions E and hence rest-masses, while main- 
taining the same radial extension. The result contained in equation 
<B4t states, therefore, that it is possible to build tori with largely 
different rest-mass densities and yet have them all share the same 
sound velocity. 

Two additional remarks are worth making. The first one is that 
while this result has been derived for a disc with a nonzero verti- 
cal extension, it applies also for a vertically integrated model. The 
second remark is that while we have been here concentrating on 
relativistic models, the same result can be shown to apply in New- 
tonian physics when the central black hole is replaced by a generic 
spherically symmetric gravitational potential. 



APPENDIX B: ON THE SOUND SPEED IN ROTATING 
POLYTROPIC TORI 

We here show that a non-self gravitating, polytropic fluid configu- 
ration orbiting in hydrostatic equilibrium around a Schwarzschild 
black hole has a sound velocity which is invariant under changes 
of the polytropic constant k. To prove this consider the equation of 
hydrostatic equilibrium <28> which, in the simplified metric <22t . 
can be rewritten as 



(Bl) 
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